CTG data analysis of HC1092 human organoid screen
24 October 2017
1 Introduction
12 organoid lines were screened with about 70 compounds (5 concentrations) of the clinical cancer panel. After 7 days total (4 days of treatment) the organoids we lysed and a ctg assay was performed. The experiment was conducted in 2 replciates.
1.1 Prepare data analysis
Load necessary packages for data analysis
#library(tidyverse)
library(readr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(corrr)
library(pheatmap)
library(gridExtra)
library(EBImage) #cite
library(HTSvis) #cite
library(PharmacoGx) #cite
library(growthcurve) #cite
library(jsonlite)
library(httr)
library(openxlsx)
library(preprocessCore)
library(stringr)
library(GRmetrics)
library(scales)Load clinical data. The primary location of tumors is divided into three sights.
#import clinical data
cohort <- read_excel("~/promise/local_data/clinical_data/Clin_Data_Basic-cohort.xlsx") %>%
separate(ID, c("p", "no.line"), sep = 1) %>%
mutate(line = paste0("D", no.line, "T01")) %>% #only works if there were no re-biopsies
select(-p, -no.line)## Warning in strptime(x, format, tz = tz): unknown timezone 'default/Europe/
## Berlin'cohort_short <- cohort %>%
#select(line, Location) %>%
select(line, age, Location, `T`, N, M, Stage, G) %>%
mutate(Patient = line) %>%
select(-line)
cohort_short$Location[cohort_short$Location == "Asc"] <- "Right"
cohort_short$Location[cohort_short$Location == "Sigm"] <- "Left"
cohort_short <- cohort_short %>%
mutate(Location = factor(Location, levels = c("Right", "Left", "Rectum")))1.2 Normalize CTG photon counts by a plates DMSO controls
1.2.1 Validation of different normalization strategies
Appropriate normalization of CTG photon counts is key to ensure robust quantification of drug response. Many different strategies are available to remove potential technical bias from the data and some of them might be more appropriate than others depending on the kinds of bias present in the data. Hence, we examine carefully the kinds of technical bias present in our data set and employ different normalization strategies in order to correct these. We assess and compare performance of these different strategies by looking at various controls.
1.2.1.1 Assessment of technical bias present in the data
##read ctg data files
ctg_data <- lapply(list.files('~/promise/rsync_data/ctg_data/HC1092/', full.names=T),
function(f) read_tsv(f, col_names=F, col_types = 'cci') %>%
`colnames<-`(c('screen', 'well', 'pcount'))) %>% bind_rows()For now we load again the plate annotation and merge that to the data to know which drug at which concentration went where.
## read excel file, convert to tibble
plate_anno <- tbl_df(read.xlsx('~/promise/rsync_data/layouts/raw_layout_files/Clin_Cancer_Panel_V170511.xlsx'))
## rename drug and well columns
colnames(plate_anno)[c(1, 6)] <- c('drug', 'well')
## join with ctg data
ctg_data <- ctg_data %>% inner_join(plate_anno) %>%
extract(well, c('row', 'col'), regex='([A-P])(\\d{2})', remove=F) %>%
mutate(row=match(row, LETTERS), col=as.integer(col))There was a technical issue with one of the plates. Wells have to be flagged accordingly.
ctg_data <- ctg_data %>%
mutate(pcount = ifelse(screen == '170502_NR_M2_D004T01P006L08' &
(row %in% seq(1, 15, by=2)) &
(col %in% 1:13), NA, pcount))Now we can get an idea of what the data looks like. First we look at only one plate at a time and look at the photon counts of DMSO controls to see how they vary and if there might be a positional bias that needs correction.
## box plot of DMSO photon counts
ctg_data %>% filter(drug == 'DMSO') %>%
ggplot(aes(rack.concentration, pcount)) + geom_jitter(width=0.2) +
facet_wrap(~screen) +
geom_boxplot(alpha=0) + theme_classic()
## heat map of DMSO photon counts, ignoring other values
p <- ctg_data %>% filter(screen == unique(ctg_data$screen)[3])
p %>% mutate(pcount = ifelse(drug != 'DMSO', -1000, pcount)) %>%
dplyr::select(row, col, pcount) %>% spread(col, pcount) %>%
dplyr::select(-row) %>% as.matrix() %>%
pheatmap(cluster_rows=F, cluster_cols=F)
DMSO controls are scattered across the plate and there is no obvious spacial bias with the exception, maybe, of the control at the bottom right of the plate. However, some of the controls clearly failed and have to be flagged form the analysis. What is concerning is that these controls are all in the same area. Maybe plotting the plate as a whole reveals a spacial bias across all drugs.
p %>% dplyr::select(row, col, pcount) %>% spread(col, pcount) %>%
dplyr::select(-row) %>% as.matrix() %>%
pheatmap(cluster_rows=F, cluster_cols=F)
I would say that there is definitely a spacial bias where photon counts get higher towards the bottom and the maybe edges of the plate.
## scatter plot with robust loess fit for rows.
ctg_data %>% ggplot(aes(row, pcount)) + geom_point() +
geom_smooth(method.args=list(family='symmetric')) +
facet_wrap(~screen) + theme_classic()
## scatter plot with robust loess fit for columns
ctg_data %>% ggplot(aes(col, pcount)) + geom_point() +
geom_smooth(method.args=list(family='symmetric')) +
facet_wrap(~screen) + theme_classic()
1.2.1.2 Normlization of photon count data
Looking at above plots this seems to be true for rowwise effect but seems to be negligible on the x axis of the plates. Hence, in order we correct for row-wise spacial bias we apply a loess-normalization to the data.
## split data by plate and apply loess normalization
ctg_loess <- ctg_data %>% split(.$screen) %>% lapply(function(s){
## loess fit. family is 'symmetric' to be robust to outliers
fit <- loess(pcount ~ row, data=s, family='symmetric')
## apply normalization
s %>% mutate(norm_fac = predict(fit, row),
pcount_norm = pcount - (norm_fac - median(norm_fac)))
}) %>% bind_rows()## plot row vs photon count of normalized values
ctg_loess %>% ggplot(aes(row, pcount_norm)) + geom_point() +
geom_smooth(method.args = list(family='symmetric')) + theme_classic() +
facet_wrap(~screen)
Next we can investigate plate specific biases and if correcting for these would actually change anything compared to just normalizing to the DMSO controls on every plate individually. We compare the effects by looking at the drug Bortezomib which we expect to be toxic regardless of the organoid line.
Of note, the 32 DMSO controls yield a reference median which is, across all plates, more robust and greater in relative size than a plate-wise median. Hence the DMSO control dependent median should be used for calibration of plates.
## no cross-plate normalization
ctg_loess %>%
left_join(ctg_loess %>% filter(drug == 'DMSO') %>% group_by(screen) %>%
summarise(ctrl = median(pcount_norm, na.rm=T)) %>% ungroup()) %>%
group_by(screen) %>% mutate(rv = pcount_norm/ctrl) %>% ungroup() %>%
filter(drug %in% c('DMSO', 'Bortezomib')) %>%
ggplot(aes(screen, rv, colour=drug)) + geom_point() +
facet_wrap(~rack.concentration) + theme_classic() +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
Without cross-plate normalization we can see that the drug is toxic, especially at higher concentrations, but we can definitely see some variation. I’m not sure at this point whether this could be a real effect or whetehr there is variation for different reasons (such as varying growth rate) that is better normalized away. Let’s see if a cross-plate normalization changes anything.
## median photon count across all dmso controls of all plates
all_dmso <- ctg_loess %>% filter(drug == 'DMSO') %>% .$pcount_norm %>% median(na.rm=T)
## cross-plate normalization
ctg_pn <- ctg_loess %>%
left_join(ctg_loess %>% filter(drug == 'DMSO') %>% group_by(screen) %>%
summarise(ctrl = median(pcount_norm, na.rm=T)) %>% ungroup()) %>%
group_by(screen) %>% mutate(pcn = pcount_norm * (all_dmso/ctrl)) %>% ungroup()
## plot bortezomib
ctg_pn %>% dplyr::select(-ctrl) %>%
left_join(ctg_pn %>% filter(drug == 'DMSO') %>% group_by(screen) %>%
summarise(ctrl = median(pcn, na.rm=T)) %>% ungroup()) %>%
group_by(screen) %>% mutate(rv = pcn/ctrl) %>% ungroup() %>%
filter(drug %in% c('Bortezomib', 'DMSO')) %>%
ggplot(aes(screen, rv, colour=drug)) + geom_point() +
facet_wrap(~rack.concentration) + theme_classic() +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
So it seems that a cross plate normalization doesn’t actually make a difference, indicating that the growth rates might be affecting the variation in phenotype mostly. Let’s see if this could be adjusted using a quantile normalization.
## generate a matrix with a column per sample and a row per well
ctg_mat <- ctg_loess %>%
left_join(ctg_loess %>% filter(drug == 'DMSO') %>% group_by(screen) %>%
summarise(ctrl = median(pcount_norm, na.rm=T)) %>% ungroup()) %>%
group_by(screen) %>% mutate(pcount = pcount / ctrl) %>% ungroup() %>%
dplyr::select(screen, well, pcount) %>%
spread(screen, pcount) %>% data.frame() %>% `rownames<-`(.$well) %>%
dplyr::select(-well)
## boxplot
boxplot(ctg_mat)
## remember rownames and colnames
rn <- rownames(ctg_mat); cn <- colnames(ctg_mat)
## quantile normalize
qnorm_mat <- ctg_mat %>% as.matrix() %>% normalize.quantiles()
rownames(qnorm_mat) <- rn; colnames(qnorm_mat) <- cn
## boxplot of normalized matrix
boxplot(qnorm_mat)
## bring back together with the meta data
ctg_qnorm <- ctg_loess %>% dplyr::select(-c(pcount, pcount_norm, norm_fac)) %>%
inner_join(qnorm_mat %>% data.frame() %>% mutate(well = rownames(.)) %>%
tbl_df %>% gather(screen, pc_norm, -well) %>%
mutate(screen = gsub('^X', '', screen)))
## look at bortezomib
ctg_qnorm %>% filter(drug %in% c('Bortezomib', 'DMSO', 'Gefitinib')) %>%
ggplot(aes(screen, pc_norm, colour = drug)) + geom_point() +
facet_wrap(~rack.concentration) + theme_classic() +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
This looks reasonable to me. As a final check we can look at the 5-10 overall most toxic drugs and checkif there are differences between the lines (expecting that there are none).
## select 5 most toxic drugs
dtox <- ctg_qnorm %>% group_by(drug) %>%
summarise(medpc = median(pc_norm, na.rm=T)) %>% ungroup() %>%
arrange(medpc) %>% .$drug %>% .[1:5]
## boxplot of these drugs across lines and concentrations
ctg_qnorm %>% filter(drug %in% dtox) %>%
ggplot(aes(screen, pc_norm)) + geom_boxplot() +
facet_wrap(~rack.concentration) + theme_classic() +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
I do not see anything that seems to really be strikingly off so I think it’s ok to precede and calculate IC50/AUC values. Viability values of > 100 % do not really make sense as we wouldn’t expect drug treatment to result in a gain of viability phenotype. Hence, we clip the values at 1 to achieve better AUC/IC50 models.
ctg_qnorm <- ctg_qnorm %>% mutate(pc_norm = ifelse(pc_norm > 1, 1, pc_norm))
## histogram of normalized viability
hist(ctg_qnorm %>% .$pc_norm)
1.2.2 Quantification of drug response
Different ways of quantifying drug response might yield different response. We first use the PharmacoGx package to calculate AUCs based on both the actual data and the fitted curve (I would assume the first to be more reasonable).
## takes a while
auc_pgx <- ctg_qnorm %>%
## separate annotation of screen name
separate(screen, c("date", "operator", "mithras", "full_barcode"), remove = FALSE) %>%
separate(full_barcode, c("line", "plate", "library"), remove = FALSE, sep = c(7, 11)) %>%
group_by(screen, drug) %>% arrange(rack.concentration) %>%
## calculate auc both for fitted curve and actual values
mutate(s_AUC_fit = computeAUC(rack.concentration, pc_norm, verbose=F,
viability_as_pct = F, area.type = 'Fitted'),
s_AUC_actual = computeAUC(rack.concentration, pc_norm, verbose=F,
viability_as_pct = F, area.type='Actual')) %>%
group_by(line, drug) %>% arrange(rack.concentration) %>%
## calculate auc both for fitted curve and actual values
mutate(l_AUC_fit = computeAUC(rack.concentration, pc_norm, verbose=F,
viability_as_pct = F, area.type = 'Fitted'),
l_AUC_actual = computeAUC(rack.concentration, pc_norm, verbose=F,
viability_as_pct = F, area.type='Actual')) %>%
ungroup()
save(auc_pgx, file = "auc_pgx.Rdata")We plot a scatter plot to see if both ways of calculating the AUC agree.
load("~/promise/auc_pgx.Rdata")
auc_pgx %>% ggplot(aes(s_AUC_fit, s_AUC_actual)) + geom_point() + theme_classic()
auc_pgx %>% ggplot(aes(l_AUC_fit, l_AUC_actual)) + geom_point() + theme_classic()It seems that they agree very well, with the exception of some outliers generated by cases with NAs. Looking through the results manually it seems that the values estimated by the fitted curve are more reasonable. It would be interesting to compare these results with the restults of a different AUC quantification, for example using a multilevel model.
1.2.3 Evaluate growth-rate adjusted compound activity
Next to quantile normalization, two recent papers suggest growth rate adjusted auc calculation to account for differences in doubling times of cell lines. First I load growth rate data from separate experiments and perform first annotation of the dataset. I filter the proliferation dataset further by removing empty wells.
#define datasets to load ctg data into R
filelist <- list.files("~/promise/rsync_data/ctg_data/proliferation_true", pattern = '_Prolif_', recursive = TRUE, full.names = TRUE)
#define a function to load proliferation ctg files into R once they match the given definitions
load_delim <- function(full.name){
read_delim(full.name,
"\t", escape_double = FALSE, col_names = FALSE,
trim_ws = TRUE) %>%
mutate(plate_name = full.name %>% str_split(., "/") %>% .[[1]] %>% .[length(.)] %>% str_sub(.,1,-5))
}
#load ctg data
ctg_prolif_raw <- lapply(filelist, load_delim) %>% bind_rows() %>% `colnames<-` (c('original_name', 'well_id_384', 'photons', "plate_barcode")) %>%
#ugly processing -it works
separate(plate_barcode, c("tmp.plate_barcode", "tmp.mu"), sep = -7, remove = FALSE) %>%
mutate(tmp.mu = substr(tmp.mu, 2,6)) %>%
separate(tmp.mu, c("mithras", "user"), sep = "_") %>%
separate(tmp.plate_barcode, c("date", "experiment", "timepoint", "lines" ), sep= "_", extra = "merge") %>%
mutate(no_lines = str_count(lines, "D0")) %>%
separate(well_id_384, c("letter", "number"), sep = 1, remove = FALSE) %>%
#data from seeding timepoints was generated in 96 well plates. For now it will be removed from the dataset to enable standardized analysis
filter(! timepoint %in% c("seeding")) %>%
#on three dates (170613, 170724 and 171010) not the complete plate was filled with organoids. Here the unseeded wells are removed from the dataset.
filter(!(date %in% c("170613") & number %in% 13:24))%>%
filter(!(date %in% c("170724") & number %in% 21:24))%>%
filter(!(date %in% c("171010") & number %in% 19:24))%>%
#one single well in the dataset has a photon count of exactly 0 (second lowest value being 658). It is set to NA
mutate(photons = replace(photons, which(photons == 0), NA))
#the data has the following distribution
ctg_prolif_raw %>%
ggplot(aes(photons)) +
geom_density(adjust = 0.3) +
scale_x_continuous(limits = c(0,50000)) +
theme_minimal() +
ggtitle("Distribution of Photon counts for all Proliferation Experiments")
I complete the annotation of the dataset by linking organoid lines to photon counts and time.
Gain a first overview about the proliferation dataset. Proliferation data for lines D015T01 and D030T01 appear discordant. Proliferation data for D019T01, D021T01 and D022T01 appear to follow a linear growth pattern
ctg_prolif %>%
group_by(date, timepoint_num, line) %>%
summarise(median = median(photons, na.rm = TRUE),
mad = mad(photons, na.rm = TRUE)) %>%
ggplot(aes(timepoint_num, median, color = date)) +
geom_point() +
geom_path(aes(group = date)) +
facet_wrap(~line) +
theme_minimal() +
ggtitle("Proliferation curves for 12 organoid lines") +
ylab("Median photon count for each timepoint") +
xlab("Time (d)")
I like to estimate the growth rate of each line in the initial phase of the prolifereation experiment. During this phase all lines show a reproducible growth rate after photon counts have been log2 transformed.
ctg_prolif %>%
select(timepoint_num, date, line, photons, photons_log2) %>%
group_by(date, timepoint_num, line) %>%
ggplot(aes(timepoint_num, photons_log2, group = date)) +
geom_jitter(alpha = 0.2, width = 0.2) +
stat_growthcurve(model = "linear") +
facet_wrap(~line) +
theme_minimal() +
scale_x_continuous(limits = c(2,4), breaks = c(2,4)) +
ggtitle("Proliferation curves for 12 organoid lines over 4 days") +
ylab("Log2 of median photon count for each timepoint") +
xlab("Time (d)") +
geom_smooth(alpha = 0.5)## Warning: Removed 1358 rows containing non-finite values (stat_growth_curve).## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## Warning: Removed 1358 rows containing non-finite values (stat_smooth).## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.6674e-31## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.6674e-31## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.6674e-31## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.6674e-31## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.6674e-31## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.6674e-31## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.99## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.6674e-31## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 1.99## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.6674e-31## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401## Warning: Removed 2926 rows containing missing values (geom_point).
Now I used the transformed photon counts to estimate a doubling rate by fitting linear models to each dataset (line and replicate).
#if multiple timepoints are estimated only these cell lines can be used, since they show log2 linear growth over all measured timepoints
picks <- c("D004T01", "D019T01", "D021T01", "D022T01")
fit <- ctg_prolif %>%
select(timepoint_num, date, line, photons_log2) %>%
#filter(timepoint_num %in% c(2,4)) %>% #there is no difference when starting the model at t=2/t=0
group_by(line, date) %>%
do(fit_growth(df = ., time = timepoint_num, data = photons_log2, model = "linear") %>% broom::tidy()) %>%
group_by(line, date) %>%
summarise(`72` = estimate[1]+estimate[2]*3,
#`96` = estimate[1]+estimate[2]*4,
#`168` = estimate[1]+estimate[2]*7,
m = estimate[2],
b = estimate[1]) %>%
gather(time, photons_log2, -m, -b, -line, -date) %>%
mutate(time = time %>% as.numeric(),
photons = 2^photons_log2,
drug = "DMSO",
full_barcode = paste0("PROLIF", line, date)) I use the fitted data to estimate the treatment date (day 3) of the screen. With this data at hand we can calculate growth rate corrected AOC and IC50 values. Moreover, I tested that the estimation of doubling time is not considerably improved if positional effects are respected and every well is handled as a separate experiment.
ctg_gr <- ctg_loess %>%
mutate(time = 168) %>%
dplyr::rename(concentration = rack.concentration,
photons = pcount) %>%
## separate annotation of screen name
separate(screen, c("date", "operator", "mithras", "full_barcode"), remove = FALSE) %>%
separate(full_barcode, c("line", "plate", "library"), remove = FALSE, sep = c(7, 11)) %>%
select(full_barcode, line, drug, time, concentration, photons) %>%
full_join(fit %>% select(line, time, photons, drug, full_barcode)) %>%
#only lines which were still in a log phase of growth on the measured timepoints can be used for estimating the needed timepoints t=3 and t=7
#filter(line %in% picks) %>%
filter(line != "D015T01") %>%
dplyr::rename(cell_count = photons,
cell_line = line,
treatment = drug) %>%
as.tibble() %>%
arrange(time)
# ctg_gr %>%
# filter(cell_line %in% picks & treatment == "DMSO" & time == 168) %>%
# ggplot(aes(full_barcode, cell_count)) +
# geom_boxplot() +
# facet_grid(~cell_line) +
# theme_minimal() +
# theme(axis.title.x=element_blank(),
# axis.text.x=element_blank(),
# axis.ticks.x=element_blank()) +
# ggtitle("Screen replicates and estimated viability at day 7") +
# ylab("Photon counts for DMSO treated wells")
drc_output <- ctg_gr %>%
select(-full_barcode) %>%
mutate(concentration = if_else(treatment == "DMSO", "0", concentration),
treatment = if_else(treatment == "DMSO", "-", treatment),
concentration = as.numeric(concentration)) %>%
mutate(time = time - 72) %>%
GRfit(., case = "C", groupingVariables = c('cell_line','treatment'))
GRdrawDRC(drc_output, experiments = c(paste0(ctg_gr$cell_line %>% unique(), " ", c("D010T01 Gefitinib", "D027T01 Gefitinib"))))
GRgetMetrics(drc_output) %>%
filter(pval_GR < 1) %>%
select(cell_line, treatment, GR_AOC) %>%
as.tibble() %>%
ungroup() %>%
dplyr::rename(line = cell_line,
drug = treatment,
aoc.mean = GR_AOC) %>%
filter(!drug %in% c( "DMSO")) %>%
#filter(!line %in% c("D015T01", "D010T01")) %>%
#filter(!grepl(.$drug, pattern = 'mit')) %>%
spread(line, aoc.mean) %>%
na.omit() %>%
as.data.frame() %>%
remove_rownames() %>%
column_to_rownames('drug') %>%
pheatmap(
#scale = "row",
cluster_rows = TRUE,
show_rownames = TRUE,
show_colnames = TRUE,
color = colorRampPalette(c('#ca0020','#f4a582','#f7f7f7', '#92c5de','#0571b0'))(5),
#color = colorRampPalette(c('#f9f902','#000000','#3701f9'))(10),
#color = c('#ca0020','#f4a582','#f7f7f7','#92c5de','#0571b0'),
annotation_col = col,
annotation_colors = ann_colors,
#annotation_row = row
#cutree_cols = 2,
cutree_rows = 4,
na_col = "grey"
) Compare AUC by the GRmetric package with the PharmacoGx analysis by correlating responses. First fit AUC data irrespective of timepoints.
compare_methods <- GRgetMetrics(drc_output) %>%
as.tibble() %>%
ungroup() %>%
dplyr::rename(line = cell_line,
drug = treatment,
aoc.mean = GR_AOC) %>%
select(line, drug, aoc.mean) %>%
mutate(aoc.mean.scale = )
left_join(auc_pgx %>%
select(line, drug, l_AUC_fit, l_AUC_actual, s_AUC_fit, s_AUC_actual) %>%
distinct(line, drug, .keep_all = TRUE))
compare_methods %>%
ggplot(aes(l_AUC_fit, l_AUC_actual)) + geom_point() + theme_classic()
compare_methods %>%
ggplot(aes(s_AUC_fit, l_AUC_fit)) + geom_point() + theme_classic()
compare_methods %>%
ggplot(aes(aoc.mean, l_AUC_fit)) + geom_point() + theme_classic()1.2.4 Illustrate differences in growth rates
For future analysis it might be interesting to know the growth rate of each line as an additional feature. I will calculate a growth rate and attach it to the annotation data which we have previously defined for each organoid line.
fit %>%
group_by(line) %>%
mutate(m_mean = mean(m)) %>%
arrange((m_mean)) %>%
ungroup(line) %>%
mutate(line = factor(line, levels = line %>% unique())) %>%
ggplot() +
geom_point(aes(x=line, y=m_mean), col="black", size=5, alpha = 0.9) + #"tomato2"
geom_point(aes(x=line, y=m), col="grey", size=2, alpha = 0.9) + # Draw points
geom_segment(aes(x=line,
xend=line,
y=min(m),
yend=max(m)),
linetype="dashed",
size=0.1) + # Draw dashed lines
labs(title="Organoid Growth Rate",
subtitle="Cell doublings inbetween 48 and 96 hours after seeding") + #caption
coord_flip() +
theme_classic() +
xlab("Organoid line") +
ylab("Cell doublings in 48h")
compound_of_interest <- c("Trifluoridin/Tipiracil")
dodge <- position_dodge(width=0.15)
ctg_qnorm %>%
separate(screen, c("date", "operator", "mithras", "full_barcode"), remove = FALSE) %>%
separate(full_barcode, c("line", "plate", "library"), remove = FALSE, sep = c(7, 11)) %>%
filter(drug %in% compound_of_interest) %>%
mutate(line = as.factor(line),
drug = factor(drug, levels = compound_of_interest),
concentration_factor = as.numeric(rack.concentration)) %>%
group_by(line, drug, concentration_factor) %>%
dplyr::summarise(mean_photons = mean(pc_norm, na.rm = TRUE),
sd_photons = sd(pc_norm, na.rm = TRUE),
range_low_photons = range(pc_norm)[1],
range_high_photons = range(pc_norm)[2]) %>%
ggplot(aes(y = mean_photons, x = concentration_factor)) +
geom_point(position = dodge, stat = "identity", aes( colour = line), size = 2) +
geom_path(aes( colour = line), alpha = 0.6, position = dodge, size = 2) +
#geom_errorbar(aes(ymin=mean_photons-sd_photons, ymax=mean_photons+sd_photons, group = line), width=0, position = dodge, size = 1, alpha = 0.5) +
#geom_ribbon( aes(x=concentration_factor, y=mean_photons, ymin=mean_photons-sd_photons, ymax=mean_photons+sd_photons), alpha=0.2) +
scale_x_log10(breaks = c(7,14)) +
ylab("Photon count normalized to control") +
facet_wrap(~drug) +
xlab("Concentration factor 5-fold, mean of n=2") +
scale_colour_brewer(palette = "Set3") +
theme_minimal()2 Session info
sessionInfo()## R version 3.4.0 (2017-04-21)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.13
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 tools stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2 GRmetrics_1.2.0
## [3] SummarizedExperiment_1.6.5 DelayedArray_0.2.7
## [5] matrixStats_0.52.2 Biobase_2.36.2
## [7] GenomicRanges_1.28.6 GenomeInfoDb_1.12.3
## [9] IRanges_2.10.5 S4Vectors_0.14.7
## [11] BiocGenerics_0.22.1 preprocessCore_1.38.1
## [13] openxlsx_4.0.17 httr_1.3.1
## [15] jsonlite_1.5 growthcurve_1.0.0.9000
## [17] PharmacoGx_1.6.1 HTSvis_1.0.2
## [19] readxl_1.0.0 shiny_1.0.5
## [21] gtools_3.5.0 DT_0.2
## [23] gplots_3.0.1 scales_0.5.0.9000
## [25] RColorBrewer_1.1-2 ggvis_0.4.3
## [27] reshape2_1.4.2 shinyjs_0.9.1
## [29] data.table_1.10.4-2 stringr_1.2.0
## [31] tibble_1.3.4 EBImage_4.18.3
## [33] gridExtra_2.3 pheatmap_1.0.8
## [35] corrr_0.2.1 ggplot2_2.2.1.9000
## [37] tidyr_0.7.2 dplyr_0.7.4
## [39] readr_1.1.1 BiocStyle_2.4.1
##
## loaded via a namespace (and not attached):
## [1] NISTunits_1.0.1 backports_1.1.1 fastmatch_1.1-0
## [4] drc_3.0-1 sm_2.2-5.4 plyr_1.8.4
## [7] igraph_1.1.2 lazyeval_0.2.0 splines_3.4.0
## [10] BiocParallel_1.10.1 SnowballC_0.5.1 TH.data_1.0-8
## [13] digest_0.6.12 htmltools_0.3.6 tiff_0.1-5
## [16] gdata_2.18.0 magrittr_1.5 cluster_2.0.6
## [19] limma_3.32.10 celestial_1.4.1 sandwich_2.4-0
## [22] piano_1.16.4 jpeg_0.1-8 colorspace_1.3-2
## [25] RCurl_1.95-4.8 lme4_1.1-14 bindr_0.1
## [28] survival_2.41-3 zoo_1.8-0 glue_1.1.1
## [31] gtable_0.2.0 zlibbioc_1.22.0 XVector_0.16.0
## [34] MatrixModels_0.4-1 car_2.1-5 maps_3.2.0
## [37] abind_1.4-5 SparseM_1.77 mvtnorm_1.0-6
## [40] relations_0.6-7 miniUI_0.1.1 Rcpp_0.12.13
## [43] plotrix_3.6-6 viridisLite_0.2.0 xtable_1.8-2
## [46] foreign_0.8-69 mapproj_1.2-5 htmlwidgets_0.9
## [49] fgsea_1.2.1 pkgconfig_2.0.1 nnet_7.3-12
## [52] locfit_1.5-9.1 tidyselect_0.2.2 labeling_0.3
## [55] rlang_0.1.2 munsell_0.4.3 cellranger_1.1.0
## [58] downloader_0.4 broom_0.4.2 evaluate_0.10.1
## [61] fftwtools_0.9-8 yaml_2.1.14 knitr_1.17
## [64] caTools_1.17.1 purrr_0.2.4 RANN_2.5.1
## [67] nlme_3.1-131 mime_0.5 quantreg_5.33
## [70] slam_0.1-40 pracma_2.0.7 compiler_3.4.0
## [73] pbkrtest_0.4-7 plotly_4.7.1 png_0.1-7
## [76] marray_1.54.0 stringi_1.1.5 lattice_0.20-35
## [79] Matrix_1.2-11 psych_1.7.8 nloptr_1.0.4
## [82] magicaxis_2.0.3 bitops_1.0-6 httpuv_1.3.5
## [85] R6_2.2.2 bookdown_0.5 KernSmooth_2.23-15
## [88] lsa_0.73.1 codetools_0.2-15 MASS_7.3-47
## [91] assertthat_0.2.0 rprojroot_1.2 withr_2.0.0
## [94] mnormt_1.5-5 multcomp_1.4-7 GenomeInfoDbData_0.99.0
## [97] mgcv_1.8-22 hms_0.3 quadprog_1.5-5
## [100] grid_3.4.0 minqa_1.2.4 rmarkdown_1.6
## [103] sets_1.0-17